#ifndef _SYSTEM_M_H_
#define _SYSTEM_M_H_
#include<string>
#include"help.h"
#include"base.h"
#include"orbital.h"
#include"twobodystate.h"
#include"system_ho.h"

using std::string;


  ///orbitals with par,mm,tz
class Orbital_M:public Orbital_parmmtz
{
 public:
  typedef double DataType;  
  Orbital_M(){}
 Orbital_M(int _par,int _mm,int _tz,DataType _e):Orbital_parmmtz(_par,_mm,_tz),e(_e){}

  DataType e;
};

  
  
  ///system hold interaction in M scheme
class System_M:public System<Orbital_M,TwoBodyState_TzParM,double>
{
 public:
  typedef double DataType;
  typedef System_Tz_HO VSystem;//interaction system
  typedef Orbital_M OrbitalType;
  typedef TwoBodyState_TzParM TwoBodyStateType;
  typedef TwoBodyMatElement<DataType> TwoBodyMatEle;
  const static int MatEledim;

 System_M(VSystem *Ptr):pSystem(Ptr){}
      
  void setup()
  {
    setupOrbitals();
    setupTwoBodyStates();
    //setupTwoBodyMat();
  }
  void setupOrbitals();
  void setupTwoBodyStates();
  void setupTwoBodyMat();

  void setupIndiceIn2BStates();
  vector<int> IndiceIn2BStates;

  void get2BmatAt(int a,int b,int c,int d,TwoBodyMatEle &MatEle) const;
  void get2BmatOnNeed(int a,int b,int c,int d,TwoBodyMatEle &MatEle) const;
      
  DataType get2BV(int A,int a,int b,int c,int d) const
  {
    TwoBodyMatEle MatEle;
    //get2BmatAt(a,b,c,d,MatEle);
    get2BmatOnNeed(a,b,c,d,MatEle);
    if(MatEle.empty()) return 0;
    return MatEle[0]-MatEle[3]*pSystem->hbar_omega/A;
  }

  DataType get1B(int A,int a,int b) const
  {
    if(Orbitals[a].isInAGroupWith(Orbitals[b]))
      return pSystem->get1B(A,Orders[a],Orders[b]);
    else
      return 0;
  }
      
  double jplus(int a,int b) const;///< J+ matrix element between two orbital a and b
  double jminus(int a,int b) const;///< J+ matrix element between two orbital a and b
      
  void printOrbitals(int num=-1) const;
  void printTwoBodyStates() const;
      
  VSystem * pSystem;
  vector<int> Orders;//split from which orbital in pSystem->Orbitals, the main shell of each orbital

 private:
  int mapab(int a,int b) const//a<b,a map from a<b matrix to one dimension vector
  {
    int totalOrbitals=Orbitals.size();
    return a*totalOrbitals+b-(a*(a+3))/2;
  }
};

#endif
